Stability of general relativistic static thick disks: the isotropic Schwarzschild thick disk 
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We study the stability of general relativistic static thick disks, as an application we consider the 
thick disk generated by applying the "displace, cut, fill and reflect" method, usually known as the 
image method, to the Schwarzschild metric in isotropic coordinates. The isotropic Schwarzschild 
thick disk obtained from this method is the simplest model to describe, in the context of General 
Relativity, real thick galaxies. The stability under a general first order perturbation of the disk is 
investigated. The first order perturbation, when applying to the conservation equations, leads to 
a set of differential equations that are, in general, not self-consistent. This opens the possibility 
of performing various kinds of perturbations to transform the resulting system of equations into 
a self-consistent system. We perform a complete classification of the perturbations as well as the 
stability analysis for all the relevant physical perturbations. We found that, in general, the isotropic 
Schwarzschild thick disk is stable under these kinds of perturbations. 
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I. INTRODUCTION 

Exact axial symmetric solutions of Einstein field equa- 
tions are a good starting point for modeling astrophysical 
applications in General Relativity. This is due to the fact 
that the natural shape of an isolated self-gravitating fluid 
is axially symmetric. In the last decades, several exact 
solutions were studied as possible galactic models. Static 
thin disk solutions were first studied by Bonnor and Sack- 
field 111 and Morgan and Morgan ^, where they con- 
sidered disks without radial pressure. Disks with radial 
pressure and with radial tension have been considered 
by Morgan and Morgan Q and Gonzalez and Letelier 
|j| , respectively. Self-similar static disks were studied by 
Lynden-BcU and Pineault Q, and Lemos €]. Moreover, 
solutions that involve superpositions of black holes with 
static disks were analyzed by Lemos and Letelier 0, H, |a| 
and Klein [l3|. Also, relativistic counter-rotating thin 
disks as sources of the Kerr type metrics were found by 
Bicak and Ledvinka [ll|. Counter-rotating models with 
radial pressure and dust disks without radial pressure 
were studied by Gonzalez and Espitia |l2 |. and Garcia 
and Gonzalez [l^l , respectively; while rotating disks with 
heat flow were studied by Gonzalez and Letelier [lj| • Fur- 
thermore, static thin disks as sources of known vacuum 
spacetimes from the Chazy-Curzon metric |l5l Il6j and 
Zipoy-Voorhees [13, [i3 metric were obtained by Bicak, 
Lynden-Bell and Katz [13. Also, Bicak, Lynden-Bell 
and Pichon 20] found an infinite number of new static 
solutions. Stationary disk models including electric fields 
[21| , niagnetic fields |23| , and both electric and magnetic 
fields ^^ have been studied. In the last years, exact 
solutions for thin disks made with single and composite 
halos of matter [23 , charged dust [23 and charged per- 
fect fluid |2g were obtained. For a survey on relativis- 
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tic gravitating disks, see [23, [23|- Most of the models 
constructed above were found using the metric to cal- 
culate its energy momentum-tensor, i.e. a inverse prob- 
lem. Several exact disk solutions were found using the 
direct method that consists in computing the metric for 
a given energy momentum tensor representing the disk 
[29; 30, .^, ^, .^, ^, .2^] . In a first approximation, the 
galaxies can be thought to be thin, this usually simplifies 
the analysis and provides very useful information. But, in 
order to model real physical galaxies the thickness of the 
disks must be considered. Exact axially symmetric rel- 
ativistic thick disks in different coordinate systems were 
studied by Gonzalez and Letelier [Sg- Also, different 
thick disks were obtained from the Schwarzschild met- 
ric in different coordinates systems with the "displace, 
cut, fill, and reflect" method [33. Recently, a relativis- 
tic generalization of the Miyamoto-Nagai potential was 
obtained i38] . 

The applicability of these disks models to any struc- 
ture found in Nature lays in its stability. The study of 
the stability, analytically or numerically, is vital to the 
acceptance of the particular model. Also, the study of 
different types of perturbations, when applied to these 
models, might give an insight on the formation of bars, 
rings or different stellar patterns. Furthermore, a per- 
turbation can cause the collapse of a stable object with 
the posterior appearance of a different kind of structure. 
An analytical treatment of the stability of disks in New- 
tonian theory can be found in [S^, [13 and references 
therein. In general, the stability of disks in General Rel- 
ativity is done in two ways. One way is to study the 
stability of the particle orbits alon g ge odesies. This kind 
of studied was made by Letelier 1411 transforming the 
Rayleigh criterion of stability [43, [3 into a general rela- 
tivistic formulation. Using this criterion, the stability of 
orbits around black holes surrounded by disks, rings and 
multipolar fields were analyzed [4j|. Also, this criterion 
was employed in 24:] to study the stability of the isotropic 
Schwarzschild thin disk, and thin disks of single and com- 
posite halos. The stability of circular orbits in stationary 



axisymmetric spacetimes are studied in |4J, |45| . More- 
over, the stability of circular orbits of the Lemos-Letelier 
solution for the superposition of a black hole and a flat 
ring are considered in |46Ll47l l4a|. Also, Bicak, Lynden- 
Bell and Katz _19] analyzed the stability of several thin 
disks without radial pressure or tension studying their ve- 
locity curves and specific angular momentum. The other 
way to study the stability of disks is to perturb its en- 
ergy momentum tensor. This way is more complete than 
the analysis of particle motions along geodesies because 
we are taking into account the collective behavior of the 
particles. However, there are few studies in the liter- 
ature performing this kind of perturbation. A general 
stability study of a relativistic fluid, with both bulk and 
dynamical viscosity, was done by Seguin [43|. He con- 
sidered the coefficients of the perturbed variables as con- 
stants, i.e. local perturbations. Usually, this condition 
is too restrictive. Stability analysis of thin disks from 
the Schwarzschild metric, the Chazy-Curzon metric and 
Zipoy-Voorhees metric perturbing their energy momen- 
tum tensor with a general first order perturbation were 
made by Ujevic and Letelier [S^l, finding that the thin 
disks without radial pressure are not stable [SJ ■ 

In order to have a general relativistic physical model 
for galaxies, we must considered, first of all, the thick- 
ness of the disk and its stability under perturbations of 
the fluid quantities. The purpose of this work is to study 
numerically the stability of the isotropic Schwarzschild 
thick disk under a general first order perturbation. The 
perturbation is done in the temporal, radial, axial and 
azimuthal components of the quantities involved in the 
energy momentum tensor of the fluid. In the general 
thick disk case, the perturbed conservation equations are 
not self-consistent because the number of unknowns is 
larger than the number of equations. In order to make 
the system of equation self-consistent we must chose a 
combination of the perturbed quantities. This opens the 
possibility to perform several types of perturbed quan- 
tities combinations that provide a self-consistent set of 
equations. In this manuscript, we study the physical ac- 
cepted perturbations. We mean by a physical pertur- 
bation the fact that a perturbation in a given direction 
of the pressure must create a perturbation in the same 
direction of the four velocity. The energy momentum 
perturbation considered in this manuscript is treated as 
"test matter", so it does not modified the background 
metric obtained from the solution of Einstein equations. 

The article is organized as follows. In Sec. II, we 
present the general perturbed conservation equations for 
the thick disk case. The energy momentum tensor is 
considered diagonal with its elements different from zero. 
Also, the energy momentum perturbation is considered 
as "test matter" . In this section we classified the possible 
perturbed quantities combinations that leads to a self- 
consistent set of the perturbed conservation equations of 
motion. In Sec. Ill, we present the thick disk considered 
for the analysis, i.e. the isotropic Schwarzschild thick 
disk. The form of its energy density and pressures, as 



well as, the restrictions that the thermodynamic quan- 
tities must obey to satisfy the strong, weak and domi- 
nant energy conditions are shown. Later, in Sec. IV, we 
perform all the physical accepted perturbations to the 
isotropic Schwarzschild thick disk to study its stability. 
Finally, in Sec. V, we summarized our results. 



II. PERTURBED EQUATIONS 

The thick disk considered is a particular case of the 
general static, axially symmetric metric 



ds' 



^a*] 



dr 



e'^^R'de' + e''^'{dR' + dz 



(1) 



where '^i, ^2 and ^3 are functions of the variables (i?, z). 
(Our conventions are: G = c = 1. Metric signature 
+2. Partial and covariant derivatives with respect to the 
coordinate x^ denoted by ,/x and ;/i, respectively.) 

In its rest frame, the energy momentum tensor of the 
fluid r^'' is diagonal with components {-p,pR,pg,Pz), 
where p is the energy density and {pR,Pe,Pz) are the 
radial, azimuthal and axial pressures or tensions, respec- 
tively. So, in this frame of reference, the energy momen- 
tum tensor can be written as 



Tf"" = pW^U" + pEXf^X" + peYf^Y" + p^.Z^'Z^ (2) 

where C/^, AT'*, F^, and Z^ are the four vectors of the 
orthonormal tetrad 



C/^-e-*ni, 0,0,0), 
e-*^ (0,1, 0,0), 



X^' 



R 



-(0,0,1,0), 



6-^^^(0,0,0,1), 



(3) 



which satisfy the orthonormal relations. Note that with 
the above deflnitions, the timelike four velocity of the 
fluid is C/^ and the quantities X^^, y , and Z^^ are the 
spacelike principal directions of the fluid. Furthermore, 
the energy momentum tensor satisfles Einstein fleld equa- 
tions, Gfj,i, = kT^v 

Due to the form of our axially symmetric metric, the 
quantities involved in the energy momentum tensor and 
the coefficients of the perturbed conservation equations 
are functions of the coordinates (i?, z) only. Let us con- 
sider a general perturbation Ap of a quantity A^ in the 
form 



A^p(t, R, e, z) = A^{R, z) + M^(t, R, 6*, z), (4) 

where A^{R^z) is the unperturbed quantity and 
(5A^(i, i?, 0, z) is the perturbation. Replacing Q for each 



quantity in the energy momentum tensor ^ and reorder- 
ing terms, we obtain that 



T^'\t, R, 0, z) = r^'^(i?, z) + 5T>"'{t, R, t 



(5) 



Hereafter, we assume that the perturbed energy mo- 
mentum tensor does not modify the background metric 
found by solving the Einstein field equation G^j/ = kT^^. 
In other words, the perturbation ST^^"^ is treated as a 
test fluid. Not considering the fluid as a test fluid is 
a second degree of approximation to the stability prob- 
lem in which the emission of gravitational radiation is 
considered. With the perturbed energy momentum ten- 
sor Tp'^ and Einstein field equations, we obtain that the 
perturbed energy momentum equations for the problem 
must obey 



Moreover, demanding the orthonormal condition for the 
perturbed tetrad we find several relations between their 
components (1411). Because of the form of the perturbed 
energy momentum equation for the axially symmetric 
metric (^, the connection coefficients are only functions 
of the coordinates {R,z). Therefore, all the coefficients 
on the perturbed conservation equations depend only on 
(R, z) and we can construct a general perturbation of the 
form 



SAf'it, R, 9, z) ^ SA^iR, z)e*('=««-™*). 



(7) 



(5T.^'' = 0. 



(6) 



These four equations can be found in Ref. 50] from Eq. 
(811) to Eq. (13II). Hereafter, the equations from Ref. 
|50| will be denoted by U. In this manuscript we are in- 
terested in linear perturbations of the energy momentum 
tensor. So, in writing the explicit form of the four equa- 
tions (|njl we discard terms of order greater or equal to 
((5^). It is explained in ^Oi|, that if we want a consis- 
tent perturbation model the tetrad must be perturbed. 



Hereafter SA^ = SA^{R, z). In this article, we are in- 
terested in perturbations of the four velocity U^ and the 
thermodynamics variables {p,pR,p0,Pz). For that rea- 
son, we see in (14II) that the t components of the tetrad 
SX*, 6Y* and (5Z* must be perturbed. Moreover, the 
tetrad perturbations 6X^ , SX^ and SY^ are not related 
whatsoever to the four-velocity or the thermodynamic 
variables perturbations and we can set them equal to 
zero. With these assumptions, the perturbation |(7J) and 
conditions (14II), the general perturbed energy momen- 
tum equations (811)- (1 311) can be cast into the particular 
form 



J 



/i = t 



fi = R. 



6U%{pU' + iiPRX"") + SU^ipU' + ^3P.^^') + '5[/^[F(i, i?, pU') + ^imPrX'' + i^¥(t, R.pnX''')] 
+5U'[ike{pU' + ^2PeY')] + <5C/^[F(t, z, pU') + Cs,.?.^" + ^^P{t. z,P.Z')] + 5p{-iwU'U') = 0, 

5pr,r{X^X^) + <5[/«[-z«;(pC/* + ^iPrX"^)] + 5p{U'U'T^) 
+SprG{R,R,X"'X'') + SpeiY^'Y^'T^g) + dp^iZ^Z^T^,) = 0, 

SU'i-wipU' + ^2PeY'^)] + SpeikeY^Y^) = 0, 

Sp,Az'z'-) + su'-Hw{pU' + ^3P,z'-)] + Spiu*u'rtt) 

+dpR{X''x''T'pp) + SpeiY'Y'Tle) + 5p,G[z, z, Z'-Z'-) = 0. 

FiI,J,K)=K,j + K{2TJj + T''^j), 
and Tp are the Christoffel symbols. 



Ai = 



p — z 



where 



(8) 

(9) 
(10) 

(11) 

(12) 
(13) 



Besides the four equations furnished by the energy mo- mentum conservation equations, TI^J" — 0, there is an- 



other important conservation equation, the equation of 
continuity, 



[nUn-,^ = 0, 



(14) 



where n is the proper number density of particles. The 
proper number density of particles n, and the total energy 
density p are related through the relation. 



p — rmib + i 



(15) 



where mb is the constant mean baryon mass and e the 
internal energy density. Multiplying Eq. (|15() by U^, 
performing the covariant derivative {; p.) and using Eq. 
p4|l . we obtain that 



(p[/^);^ = ieUn:., 



(16) 



Now, from the relation U,yT.^^ = and the energy mo- 
mentum tensor ((SJ, we obtain and expression for {pU^)-f^. 
Substituting this last expression into Eq. (|16|l we finally 
arrive at 



isU'^).^ ^ PRX''U,Xl'^+pgY''U,Y.;;^+p,Z^U,Z';'^, (17) 

which is a first order differential equation for e. There- 
fore, with e given by H17(l the equation of continuity 114|l 
is satisfied. For this reason, the continuity equation can 
be omitted for our stability analysis because, in princi- 
ple, we can always find a solution for s. Hereafter, the 
contribution of nrnf, and e to the total energy density are 
consider in p. In the case when the internal energy den- 
sity of the fluid is given, the equation of continuity must 
be considered. In other words, the knowledge of the ther- 
modynamic properties of the disk is important for actual 
applications. The thermodynamic properties of the sys- 
tem can be obtained from observations or theoretically 
from the Fokker-Planck equation, where we obtain the 
particle distribution function of the disk. Solving the 
Fokker-Planck equation is not an easy task, but some 
progress in Newtonian gravity have been done [sj . 

The four equations H8I11II contain seven independent 
unknowns, say 6U'^,SU^ ,SU^ ,6p,6pii,Spg,6pz. So, at 
this point, the number of unknowns are greater than the 
number of equations and the system is not self-consistent. 
This opens the possibility to perform different kind of 
perturbations to make the system self-consistent. In Ta- 
bleQ]we present the classification for all the possible per- 
turbations that lead to a self-consistent system of equa- 
tions. We divided the table in three sections that depend 
on the number of perturbed quantities, i.e. two, three 
and four quantities perturbations. Equation HlOf) is im- 
portant for this classification because it only involve the 
two unknowns {SU^ ,pg). Therefore, if one of them is, by 
construction, equal to zero then the other unknown is also 
zero. There is a possibility to avoid the other unknown to 



be equal zero, i.e. if the unknown coefficient is equal zero. 
This situation can occur in particular space-time points. 
We do not consider these cases. The system of equations 
for the different combinations of the perturbed quantities 
reduce to various differential equations, the form of the 
particular differential equation is written in the column 
Differential Equation. We denote by ODE# and PDE# 
an ordinary different equation of order # and a partial 
differential equation of order #, respectively. The differ- 
ential equation marked with a question mark (?) means 
that other four perturbations, that do not include the 
SU^ and Spe perturbations, are not self-consistent, but 
they might be if we assume other imposed conditions 
for the involve quantities, as for example Spfi = Spz- 
Moreover, if we assume other imposed conditions for the 
different types of perturbation shown in Tabled the dif- 
ferential equation classification may change. 

We are interested in consistent physical perturbations. 
For that reason we search for perturbations in which the 
velocity perturbation in a certain direction leads to a 
pressure perturbation in the same direction. For ex- 
ample, if we perturbed the z component of the veloc- 
ity 5U^ then we must perturb Spz to have a consistent 
physical problem. With the above criterion only four 
perturbations are allowed, these are marked with (*) in 
Table Q] These perturbations will be considered in our 
thick disk model. Furthermore, we perform the pertur- 
bation SU^, Spa, SU^, Spz with the extra imposed condi- 
tion SpR = Spz- This last perturbation are among the 
ones denoted by question mark. In this particular case, 
the system of equations reduced to a second order partial 
differential equation. 



III. ISOTROPIC SCHWARZSCHILD THICK 
DISK 

The static isotropic Schwarzschild thick disk was found 
in 37] by applying the displace, cut, fill and replace 
method to the axially symmetric metric (Q with the def- 
initions 



*i = In 



^„ = Ijn 



f 2r — m 
\2r + m J 

21nfl 



2r 



(18) 
(19) 



where r^ — R"^ + {h + 6)^, {m > 0, b > 0) are parameters, 
and ft, is a function of the coordinate z given by 



hiz) 



—z + C, z < —a, 

y4z2-|-Bz2"+2, ~a< z<a, 
z + C, z > a, 



(20) 



with 



TABLE I: Differential equation classification for all the pos- 
sible different types of perturbations. The perturbations 
marked with (*) are the ones that are physical accepted. 
Other physical accepted perturbations are marked with (?) 
but these perturbations require extra imposed conditions. 



Perturbations 


Differential Equation 




Two Perturbations 




SpR 


Spz 


PDE2 


' 


Three Perturbations 


su"" 


SpR 


5U' 


PDE2 


su"" 


SpR 


Spz 


ODEl + PDE2 


su"" 


SpR 


Sp 


ODE2 * 


w 


5W 


Spz 


PDE2 


su"" 


SW 


Sp 


PDEl 


su"" 


Spz 


Sp 


PDE2 


SpR 


SU" 


Spz 


ODEl + PDE2 


SpR 


SU' 


Sp 


PDE2 


su^ 


Spz 


Sp 


ODE2 * 




Four Perturbations 


SU" 5pe 


5[/« 


SpR 


ODE2 * 


SU" Spg 


5[/« 


3U^ 


PDEl 


SU" Spg 


JC/« 


Spz 


PDE2 


SU" Spe 


5C/« 


Sp 


ODEl 


SU" Spg 


SU" 


Spz 


ODE2 * 


SU" Spg 


SU" 


5pR 


PDE2 


SU" Spg 


SU' 


Sp 


ODEl 


SU" Spg 


SpR 


Sp 


ODEl 


SU" Spg 


5pz 


Sp 


ODEl 


Without 5C/*, 


Spg 


? 



A = 
B = 



2ri + 1 — ac 

Ana 
ac — 1 

4n(n -I- l)a2"+i ' 

a{2n + 1 + ac) 

4(n+l) ' 



where n = 1, 2, • • •, and c is the jump of the second deriva- 
tive on 2; = ±a, see J33- The energy density and the 
different pressures for this thick disk are 



64:171 



(m + 2r)' 



PR =Pe 



:[3{h + br{l-h'^) 

+r^[h'^-l + h"{h + b)]], (21) 
^^""^ -[2(/i + 6)2(1 -/,'2) 



(m + 2r)5(-m + 2r) 

+r^[h'^-l + h"{h + b)]], (22) 
64m2(l-/i'2)(/i-H6)2 



(m + 2r)5(-7Ti + 2r) 



(23) 



without losing generalization, the range of the z coordi- 
nate between —l<z<l. These means that the param- 
eter a of [23| is equal to one. To satisfy the strong energy 
condition we must have that e = p + pn + pe + Pz > 0, 
the weak energy condition requires p > and the dom- 
inant energy condition requires \pr/p\ < 1, \pe/p\ < 1 
and \pz/o'\ < 1- If we set the parameters {m,b,n,c) ex- 
plained in jSTj equal to (m = 1,6 = 2,n = l,c = 0), 
we obtain a thick disk with all energy conditions sat- 
isfied. Also, the level curves for this disk show that it 
is physically acceptable. For this set of parameters, the 
function h has the form h — i/Az"^ — 1/Sz^. We remark 
that these are not the only parameters in which the level 
curves are physically accepted. In the next section we ap- 
plied the selected perturbations of Sec. II to the isotropic 
Schwarzschild thick disk mention above and studied its 
stability. 



IV. PERTURBATIONS 

Before we start applying the different kinds of pertur- 
bations to the isotropic Schwarzschild thick disk we must 
have some considerations. Note that the exact thick disk 
metric considered is infinite in the radial direction, but 
finite in the z direction. In order to study the stability we 
need a criterion to make a finite disk, i.e. we need a cutoff 
in the radial coordinate. In Eqs. (|21|l . (|22|l and H23() . we 
see that the thermodynamic quantities decrease rapidly 
enough to define a cutoff radius. The cutoff radius Rcut is 
set by the following criterion: the energy density within 
the thick disk formed by the cutoff radius is more than 
90% of the infinite thick disk energy density. This leads 
to a cutoff radius of Rcut = 20 units. Furthermore, the 
other 10% of the energy density that is distributed from 
outside the cutoff radius to infinity can be treated, if 
necessary, as a perturbation in the outermost boundary 
condition. 



A. Perturbation with 5U , Spg, SU^, SpR 

We start perturbing the four velocity in its components 
9 and R, from the physical considerations mentioned in 
Sec. II we also expect variations in the thermodynamic 
quantities Pe andpi^;. The set of equations l|5|l- lfTT|) reduce 
to a second order ordinary differential equation for the 
perturbation Spn, say 



FaSpr^rr + FbSprr + FcSpR = 0, 



(24) 



where the prime (/) represents derivative with respect to 
the axial coordinate z. We have set in Eqs. (|21|) - (|23|l . 



where {Fa, Fb, Fq) are functions of (i?, z, w, kg), see Ap- 
pendix ^ For this particular case we have for the per- 
turbation, 5pg = —SpR. 

Note that in Eq. (|24|l the coordinate z only enters like 
a parameter. Moreover, the equation for SpR is indepen- 
dent of the parameter w, but w needs to be different from 



zero to reach that form. The second order equation H24() 
is solved numerically with two boundary conditions, one 
at _R w and the other at the cutoff radius. At i? w 
we set the perturbation Spj^ to be « 10% of the unper- 
turbed pressure pji H22f) . In the outer radius of the disk 
we set 6pr\r=r^^^ = because we want our perturbation 
to vanish when approaching the edge of the disk, and in 
that way, to be in accordance with the linear perturba- 
tion applied. We say that our perturbations are valid if 
their values are lower, or of the same order of magnitude, 
than the 10% values of the unperturbed quantities. 

In Fig. n] we present the amplitude profile of the radial 
pressure perturbation, the physical radial velocity per- 
turbation SU^ = ^/gRR^U^ and the physical azimuthal 

velocity perturbation SU^ = ^goeSU^ in the plane z — Q. 
In Fig. n we see that the perturbation 5pR decreases 
rapidly with R and has oscillatory behavior. At first 
sight, the perturbation 5pR appears to be stable for all 
R, but in order to make a complete analysis we have to 
compare at each radius the values of the different pertur- 
bations with the values of the radial pressure. For this 
purpose, we included in the same graph a profile of the 
10% value of pr. We see that the perturbations of 5pR 
for different values of kg are always lower or, at least, 
of the same order of magnitude when compared to these 
10% values. 

Our four velocity f/'' Q has only components in the 
temporal part, so we do not have values of U^ and U^ 
to make comparisons with the perturbed values 5U^ and 
5U^ of FigH For that reason we compared, in first ap- 
proximation, the amplitude profiles of these perturba- 
tions with the value of the escape velocity in the Newto- 
nian limit. In the Newtonian limit of General Relativity, 
we have the well known relation goo = ^700 + 2$. So, the 
Newtonian escape velocity Ve = •\/2|$| can be written as 



Ve 



(1 



(1 



2r> 



1/2 



(25) 



For the parameters considered, we have that at (i? — 
0, 2 = 0) and [R = 20, z = 0) the escape velocity is 
Ve — 0.8 and Ve ~ 0.31, respectively. The escape veloc- 
ity values have a small variation with z. So, with this 
criterion, we may say that the perturbations 6U^ and 
SU^ are stable because their values are always well below 
the escape velocity value. Recall that the perturbation 
SpR does not depend on the parameter w, but the per- 
turbations 6U^ and 6U^ do. In Fig. |2] we show the 
numerical solution of 5U^ and 6U^ for different values of 
the frequency w with the same wavenumber kg. We see 
that when we increase the value of the parameter w the 
perturbations become more stable. This behavior is the 
same for other values of kg . 

In this subsection we set the value of the parameter 
z = 0. We performed the same analysis for different val- 
ues of the parameter —l<z<l, and we found that 



the perturbations show the same qualitative behavior. 
Therefore, we can said that the isotropic thick disk is 
stable under perturbations of the form presented in this 
subsection. Nevertheless, if we treat the 10% of the en- 
ergy density as a perturbation in the outermost radius 
of the disk by setting Spr\r^r^^^^ — e, where e < 10% 
of PR\R=R^^t the qualitative behavior of the mode pro- 
files are the same. From these results, we can say that 
the general linear perturbation applied is stable, and for 
that reason, the perturbations do not form more complex 
structures like rings, bars or halos. Moreover, if we set 
the frequency w ^^ iw we obtain the same equation for 
the perturbation 5pR, say (|^ . In this case, the real part 
of the general perturbation lO diverge with time and the 
perturbation is not stable. These last considerations can 
be applied to every perturbation in the following subsec- 
tions. 



B. Perturbation with 5U , Spo, SU^, Spz 

In this subsection we perturb the four velocity in its 
components 9 and z, and we expect variations in the ther- 
modynamic quantities pg and pz- The set of equations 
(|5|)- (lllll reduce to a second order ordinary differential 
equation for the perturbation 5pz given by 



FaSpz 



FbSpz,z+FcSpz=0, 



(26) 



where (Fa, Fb, Fc) are functions of (i?, z,w,kg), see Ap- 
pendix^ Note that in Eq. (|26|l the coordinate R only 
enters like a parameter. Like the previous case, Eq. l|^ 
is independent of the parameter w, but in order to reach 
that form we must have w different from zero. The sec- 
ond order equation (|26l) is solved numerically with two 
boundary conditions, one in z = and the other in z = 1. 
At z = we set the perturbation 5pz to be « 10% of the 
unperturbed pressure Pz H23|l . In the outer plane of the 
disk we set 5pz\z=i = because we want our perturba- 
tion to vanish when approaching the edge of the disk, 
and in that way, to be in accordance with the linear per- 
turbation applied. 

In FigO we present the amplitude profiles of the axial 
pressure perturbation, the azimuthal pressure perturba- 
tion, the physical axial velocity 5U^ — ^/g^SU^ and the 
physical azimuthal velocity for the value of the parame- 
ter i? = 0.1. For comparison reasons, we included in the 
graph of the axial pressure perturbation the amplitude 
profile that corresponds to 10% of the value oi Pz- We 
see in this graph that the modes profiles of the axial pres- 
sure are below, or of the same order of magnitude, when 
compared to the 10% profile. In this example the modes 
with wavenumber kg — and kg = 10 have within the 
disk regions in which the profile values exceed the 10% 
profile. This perturbations are near the validity criterion 
used for the perturbations. In the axial velocity pertur- 
bations graph we see that the profiles suffer a sudden 



variation near the edge of the disk z — \. The axial ve- 
locity perturbation tends to expel the particles out from 
the disk. In our case, the particles are not allow to escape 
from the disk due to the fact that our geometrical disk is 
only valid for —1 < < 1. For that reason, the particles 
feel a kind of a wall that impedes them from going out. 
It is clear from the axial velocity perturbation profiles 
that the two nodes (w = 1, fcg = 0) and (w = l,ko = 10) 
are not stable, while the rest of the nodes are stable if 
we discard the region that is near the edge of the disk at 
z = 1. Note that the starting point of the sudden vari- 
ation of the axial velocity perturbation coincides with 
the radius in which the axial pressure perturbation ex- 
ceeds the 10% value profile of pz- Also, in Fig. |31 we 
see that the azimuthal pressure perturbations and the 
azimuthal velocity perturbations for different modes are 
stable. In the case of the azimuthal pressure, we have not 
included the 10% profile because is three order of magni- 
tude higher than the perturbation value. Note that the 
amplitude values of azimuthal velocity perturbation are 
below the escape velocity considered for the disk. 

The perturbation 5pz does not depend on the parame- 
ter w, but the perturbations 5U^ and 5U^ do depend. In 
Fig. 21 we show the profile for the mode kg = Q with dif- 
ferent values of the frequency w for the perturbation 5U^ . 
We see that when we increase the value of the parameter 
w = 1, the velocity perturbations are more stable (if we 
neglect the region near the edge of the disk at z — 1). 
This is due to the fact that the frequency w only enters 
the equation for SU'^ like a scale factor. 

We have performed the same above analysis for differ- 
ent values of the parameter R, and we found that the 
qualitative behavior is the same. With these results we 
can say, if we neglect the region near the edge of the disk 
at z = 1, that the isotropic thick disk is always stable 
for higher values of the parameter w and has some stable 
modes for lower values of w. Moreover, for lower values 
of the frequency w, some modes of our general first order 
perturbation are not stable. In these cases, more com- 
plex structures may be form, but in order to study them, 
we need a higher order perturbation formalism. 



C. Perturbation with 5U^, Spn, 5p 

In this subsection, we perturb the radial component 
of the four velocity, the radial pressure and the energy 
density of the fluid. The set of equations (|5|l- Hll|l reduce 
to a second order ordinary differential equation for the 
perturbation SpR of the form H24|) . The form of the func- 
tions {Fa, Fb,Fc) are given in Appendix lO In this case, 
the coordinate z only enters like a parameter. Due to the 
fact that we are not considering perturbations in the az- 
imuthal axis, the coefficients of the second order ordinary 
diffetential equation do not depend on the wavenumbcr 
ke ■ This second order equation is solved numerically with 
the same boundary conditions described in Sec. IV-A. 

In Fig. [Slwe present the amplitude profiles for different 



perturbation modes of the axial pressure, physical axial 
velocity and energy density for the value of the parame- 
ter z — 0. We see in the graphs, that the perturbation 
profiles decrease rapidly in few units of R. Note that the 
graphs in Fig. Oare plotted in the range [0,2], but the 
radius of the disk has been set in 20 units. We do not 
plot the 10% profile of the energy density because is at 
least two order of magnitude higher than the values de- 
picted. Also, the values of the radial perturbation modes 
profiles are well below the escape velocity. We can say 
that all these modes are highly stable. We performed the 
above analysis for different values of z and we found that 
the quantities involve have the same qualitative behav- 
ior. From these results, we can say that the general linear 
perturbation applied is highly stable, and for that reason, 
the perturbations do not form more complex structures 
like rings, bars or halos. 



D. Perturbation with SV , 5pz, Sp 

In this subsection we perturb the z component of the 
four velocity, the z component of the pressure and the 
energy density of the fluid. The set of equations l|S|l- Hll|) 
reduce to a second order ordinary differential equation 
for the perturbation 6pz of the form (j^ . The functions 
{Fa, Fb, Fc) are given in Appendix IdI Note that, like in 
Sec. IV-B, the coordinate R enters like a parameter. In 
this case, we are not considering azimuthal perturbations 
and therefore the quantities involve do not depend on 
the parameter kg. The second order equation is solved 
following the procedure of Sec. IV-B. 

In Fig. El we present the amplitude profiles of the ax- 
ial pressure perturbation, the axial velocity perturbation 
and the energy density perturbation for the value of the 
parameter R = 0.1. We see that the the pressure pertur- 
bation modes are always of the some order of magnitude 
or lower when compared to the 10% profile. Similar to 
the Sec. IV-B case, the velocity perturbation amplitude 
profiles increase abruptly near the edge of the disk at 
z = 1. For higher values of the parameter w the modes 
are stable. Even for the most stable mode there is a 
abruptly increase in the amplitude profile at the edge of 
the disk. This change in the behavior occurs very near 
z — 1 and can not be seen in the axial velocity perturba- 
tion graph. In this case we can say that the mode with 
k; = 1 is not stable. If we compare the axial velocity per- 
turbation graph of Fig. (HI with Fig. ^ we note that the 
perturbation modes considered in this subsection have 
larger regions of stability. The modes that correspond 
to the energy density perturbations are all stable. Note 
that we have not plot the 10% profile of the energy den- 
sity because is at least two order of magnitude higher 
than the values depicted. We performed the above anal- 
ysis for different values of the parameter R and we found 
that the quantities involve have the same qualitative be- 
havior. Similar to Sec. IV-B, some modes corresponding 
to lower values of the frequency w are not stable. As we 



note in Fig. (BJ this is related to the fact that our first 
order perturbation is not vaUd in some regions. In these 
cases, more complex structures may be form, but in or- 
der to study them, we need a higher order perturbation 
formalism. 



E. Perturbation with SU^, Spn, W , Spz and 

5pR = Sp:, 
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radial perturbation. The center of the disk tries to con- 
densate. So, with these considerations, we may say that 
the disk tends to form some kind of a ring around the 
center of the disk. This phenomenon appears near the 
center of the disk in a few units of R. The parameter w 
enters only like a scale factor in the differential equations 
quantities studied. 



In this subsection we perturb the radial component of 
the four velocity, the axial component of the four velocity, 
the radial pressure and the axial pressure. This kind of 
perturbation belongs to the ones marked with a question 
mark in Tabled As we said in Sec. II, we need an extra 
condition to make the set of equations self-consistent. 
In this case, we set 5pR = 6pz = Sp. Therefore, the 
set of equations l|S|)- ((TT1) reduce to a second order partial 
differential equation for the pressure perturbation Sp, say 



FaSp,rr + FbSpm + FcSp^z,+FdSp^,+FeSp = 0, (27) 

where {Fa,Fb,Fc,Fd,Fe) are functions of {R,z,w), 
see Appendix |E] The partial differential equation H27() 
is solved numerically with four boundary conditions, at 
z « —1, z « 1, i? « and R = 20. The boundary con- 
dition at z can not be set at 1 because some quantities 
diverge. As we said in Sec. IV-B, this is due to the geo- 
metrical restrictions that has the isotropic Schwarzschild 
thick disk. They are different ways in which we can 
set the boundary conditions in order to simulate various 
kinds of pressure perturbations. Here, we treat only the 
case when we have a pressure perturbation at i? « and 
along the z axis, i.e. some kind of a rod perturbation. We 
set the value of the rod pressure perturbation to be 10% 
of the axial pressure. For the rest of the boundary condi- 
tions we set the values equal to zero because we want the 
perturbation to vanish when approaching the edge of the 
disk. We chose the 10% of the value of the axial pressure 
instead of the radial pressure because it has the lowest 
value near i? w 0. In that way, the perturbation values 
are also below the 10% values of the radial pressure and 
the general linear perturbation is valid. 

In Fig. \7\ we present the perturbation amplitude for 
the pressure, the physical radial velocity and the phys- 
ical axial velocity for w = 1. We see in the pressure 
perturbation graph that the perturbation rapidly decays 
to values near zero when we move out from the center 
of the disk. In the velocity perturbations profiles we can 
see an interesting phenomenon. Note that in the lower 
domain of the disk [-1,0) the axial velocity perturbation 
is positive and in the upper domain (0,1] the axial ve- 
locity perturbation is negative. This means that due to 
the linear perturbation the disk tries to collapse to the 
plane z — 0. Now, if we look to the radial velocity per- 
turbation graph, we note that the upper and lower parts 
depart from the center of the disk due to the positive 



V. CONCLUSIONS 

In this manuscript we obtain the basic perturbation 
equations for studying the stability of thick disks. These 
equations were obtained when we applied a general first 
order perturbation to the conservation equations of mo- 
tion. The great number of unknowns make the set of 
equations not self-consistent. In order to make the sys- 
tem of equations self-consistent we must chose a combina- 
tion of the perturbed quantities. This opens the possibil- 
ity for various types of combinations. In this manuscript 
we made a classification of these possibilities and stud- 
ied the physical accepted perturbations. We can say, that 
the stability analysis we performed is more complete than 
the stability analysis of particle motion along geodesies 
because we take into account the collective behavior of 
the particles. However, this analysis can be said to be 
incomplete because the energy momentum perturbation 
tensor of the fluid is treated as a test fluid and does not 
alter the background metric. This is a second degree 
of approximation to the stability problem in which the 
emission of gravitational radiation is considered. 

The stability analyses made to the isotropic 
Schwarzschild thick disk, show that this disk is stable 
when we applied radial and azimuthal perturbations. In 
the case of axial perturbations, the nodes are not stable 
near the edge of the disk. The lack of stability is due 
to the form of the geometric thick disk considered |37j . 
i.e. it is only valid in a flxed axial interval —a<z<a, 
where a is a constant (in our case a=l). If we discard 
the region near the axial edge we can say, in general, 
that for higher values of the parameters w and kg the 
axial modes are also stable. For lower values of w, some 
axial perturbations are not stable. This is related to the 
fact that our first order perturbation is not valid in some 
regions of the considered quantities. In these cases, more 
complex structures may be form, but in order to study 
them, we need a higher order perturbation formalism. 
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APPENDIX A: COEFFICIENTS OF 
PERTURBATION TYPE A 

The general form of the functions {Fa, Fb, Fq) appear- 
ing in the second order ordinary differential equation H24|) 
are given by 

FA=Aiai, Fb ^ Aiaip. + Aia2+ A^ai, 
Fc = Aia2,R + A3a2+Aia3, (Al) 

where ai, a2 and a^, are 



where ai, a2 and aa are 



ai 



as 



Bj_ 
B2' 

£± 



B.Da - B.D, 



a2 



4^3 



B2D3 



(C2) 



where the meaning of the notation for the coefficients 
{Ai, Bi, Di) are explained in Appendix IXI 



ai 



as 



_Bj^ 
B2' 
C2D4 
C1D5' 



a2 



B5D4 — B4D5 

B^5 



(A2) 



In Eqs. IjAip and (|A2|) . we denote the coefficients of Eq. 
^ by Ai , the coefficient of Eq. (jSJ by i?i , the coefficient 
of Eq. lfTU|l by Ci, the coefficient of Eq. ((TT|l by D^, e.g., 
the first term in has the coefficient Ai multiplied by 
the factor SU^, the second term has the coefficient A2 
multiplied by the factor 6U^, etc. The explicit form of 
the above equations is obtained replacing the fluid vari- 
ables {p,PR,Pe,Pz) of the isotropic Schwarzschild thick 
disk. 



APPENDIX B: COEFFICIENTS OF 
PERTURBATION TYPE B 

The general form of the functions {Fa, Fb, Fc) appear- 
ing in the second order ordinary differential equation 126(1 
are given by 

FA=A2ai, Fb = A2ai^^+ A2a2+ A^ai, 
Fc = A2a2,z + Aias + ^5^2, (Bl) 

where ai, a2 and as are 



D2' 
C2BQ 



a2 



BqD^ — B^Dq 



(B2) 



where the meaning of the notation for the coefficients 
{Ai, Bi,Ci, Di) are explained in Appendix IXI 



APPENDIX C: COEFFICIENTS OF 
PERTURBATION TYPE C 

The general form of the functions {Fa,Fb,Fc) are 
^iven by 

FA=Aiai, Fb ^ Aiai^R + Aia2+ A^ai, 
Fc = Aia2.R + A3a2 + Aea3, (CI) 



APPENDIX D: COEFFICIENTS OF 
PERTURBATION TYPE D 

The general form of the functions {Fa,Fb,Fc) are 
given by 



FA^A2ai, Fb ^ A2ai.^ + A2a2 + A^ai, 
Fc^A2a2,,+A5a2 + Aea3, (Dl) 

where ai, a2 and as are 



ai 



as 



D2' 

El 
B3' 



BaD3 - B3D, 



a2 



3^6 



B3D2 



(D2) 



where the meaning of the notation for the coefficients 
{Ai, Bi,Di) are explained in Appendix IXI 



APPENDIX E: COEFFICIENTS OF 
PERTURBATION TYPE E 

The general form of the functions {Fa, Fb,Fc,Fd,Fe) 
appearing in the partial second order differential equation 
(|?7jl are given by 



FA = Aiai, Fb = ^iaijj-|-^ia2 -I- Asai, 
Fc ^ A2a3, Fd = A2a3^^ + A2a4 + A^as, 
Fe = Aia2,R + A2a4^z + A3a2 + A^ai, (El) 

where ai, a2, as and a^ are 



ai = 


Bl 

B2 


as = 


Dl 
D2 



012 



a^ 



B4 + Bq 

D4 + D6 

Dn ■ 



(E2) 



where the meaning of the notation for the coefficients 
{Ai, Bi, Di) are explained in Appendix IXI 
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FIG. 1: Profiles of the amplitude perturbations modes for the radial pressure, radial physical velocity and azimuthal physical 
velocity for the case when z = 0, w = 1 and ke = 0, 1, 2, 3, 4. In the pressure perturbation graph the 10% profile value of the 
radial pressure is depicted to ensure that the perturbations are in accordance with the stability criterion imposed. The value 
of the velocity perturbation profiles are well below the particles escape velocity. 
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FIG. 2: Profiles of tfie amplitude perturbations modes for the radial physical velocity and the azimuthal velocity for the case 
when z = 0, ko = 3 and w = 1, 5, 10. We note that when we increase the value of the frequency w the modes are more stable. 
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FIG. 3: Profiles of the amplitude perturbations modes for the axial pressure, azimuthal pressure, axial physical velocity and 
azimuthal physical velocity for the case when R — 0.1, w = 1 and ke = 0,10,20,30,40. In the pressure perturbation graph 
the 10% profile value of the axial pressure is depicted to ensure that the perturbations are in accordance with the stability 
criterion imposed. In the azimuthal pressure perturbation graph the 10% profile was omitted because the values are three order 
of magnitudes higher. The values of the azimuthal velocity perturbation profiles are well below the particles escape velocity 
and the modes are stable. For the axial velocity perturbation profiles we see that the geometric properties of the disk create a 
kind of wall that impedes the particles for going out. 
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FIG. 4: Profiles of the amplitude perturbations modes for the axial physical velocity for the case when 7? = 0.1, ke = and 
w = 1, 5, 10. We note that if we increase the value of the w parameter, the region in which the mode is stable increases. For 
higher values of w the modes are stable if we discard the region near z = 1. 
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FIG. 5: Profiles of the amplitude perturbations modes for the radial pressure, radial physical velocity and energy density for 
the case when a = and w = 1, 5, 10, 15, 20. As in Fig. we have depicted the 10% profile value of the radial pressure in the 
pressure perturbation graph. We see from these graphs that the pressure perturbation values are well below the 10% profile. 
In the energy density perturbation graph the 10% profile of the energy density was omitted because the values are two order of 
magnitudes higher. Also, the values of the velocity perturbations profiles are always lower than the escape velocity. Therefore, 
the modes are highly stable. 
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FIG. 6: Profiles of the ampUtude perturbations modes for the axial pressure, axial physical velocity and energy density for 
the case when R = 0.1 and w = 1,5, 10, 15, 20. To ensure that the perturbations are in accordance with the stability criterion 
imposed, we have depicted in the pressure perturbation graph the 10% profile of the axial pressure. In the energy density 
perturbation graph the 10% profile of the energy density was omitted because the values are two order of magnitudes higher. 
In the axial velocity perturbation profiles the mode with w = 1 is not stable. Note that the starting radius of the abrupt increase 
in the axial velocity perturbation amplitude coincide with the point in which the pressure perturbation becomes greater than 
the 10% profile. 
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FIG. 7: Profiles of the amplitude perturbations modes for the pressure, radial physical velocity and axial physical velocity for 
the case when w = 1. We see that the perturbation rapidly goes to values near zero when we depart from the center of the 
disk, i.e. the perturbations are stable. The parameter w enters as a scale factor in the differential equations and do not alter 
the qualitative properties of the perturbations. 



